function [pointRoot] = newtonsMethod(func, approxRoot, tol, maxIter, typeTol)
%NEWTONSMETHOD Function that determines the root of a linear function
%using Newtons Method.
%   This algorithm was adapted from Chapter 2.3 of Numerical Analysis (8th
%   edition, Burden & Faires)
%   
%   *Input Parameters*
%   func: Sym object of function (of x) to be evaluated to find root.
%   approxRoot: Approximate value of the (a) root of func.
%   tol: Tolerance limit. 
%   maxIter: Maximum number of iterations
%   typeTol: String that defines how the tolerance condition is evaluated.
%            See tolCondition function in this file for available options.
%
%   *Default options* 
%   tol = 10^-5
%   maxIter = 100
%   typeTol = 'rangeTol' -> (pointB-pointA)/2
%
%
%   ***If no root is found after maxIteration, the estimate for maxIter 
%   approximation is returned.***
%
%   ***********************************************************************
%   Author: Mathieu Boudreau, BSc, MSc, PhD Candidate (BME)
%   Institute: Montreal Neurological Institute, McGill University
%   Contact: mathieu.boudreau2 (at) mail.mcgill.ca
%   Date: July 17th 2014
%   ***********************************************************************


%% Set default input conditions, if required.
%
if nargin < 5
    typeTol = 'absMidDiff';
end

if nargin < 4
    maxIter=100;
end

if nargin < 3
    tol=10^-5;
end

%% Initialize initial conditions
%
initIter=1;

%% Run Newton's Method Algorithm
%
for iter=initIter:maxIter
    currentEstimateRoot=approxRoot - subs(func,approxRoot)/subs(diff(func),approxRoot);
    
    if tolCondition(currentEstimateRoot,approxRoot,func,typeTol)<tol
        disp(['Algorithm terminated. Tolerance condition was reached after ' num2str(iter) ' iterations.'])
        pointRoot = currentEstimateRoot; 
        return
    end
    
    approxRoot=currentEstimateRoot;
end

if iter==maxIter
    disp('Function terminated after the maximum number iterations. No solution found within tolerance.')
    pointRoot = currentEstimateRoot; 
end

end

function tolCalc=tolCondition(currentEstimateRoot,approxRoot,func,typeTol)
    %% This function chooses and evaluates the tolerance condition. 
    % One of each flavour, take your pick! 
    
    switch typeTol
        case 'absMidDiff'
            tolCalc=abs(currentEstimateRoot-approxRoot);
        case 'relMidDiff'
            tolCalc=abs(currentEstimateRoot-approxRoot)/abs(currentEstimateRoot);
        case 'funcValue'
            tolCalc=abs(subs(func,currentEstimateRoot));
    end
    
end